{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global isotropic elastic SH wave propagation using an axisymmetric approximation\n",
    "\n",
    "So far we solved the isotropic elastic equations of motion or their 2D approximations only in Cartesian coordinates. However, to describe global seismic wave propagation, we have to rewrite the governing equations in spherical coordinates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equations of motion for 3D isotropic elastic wave propagation in spherical coordinates\n",
    "\n",
    "According to [this lesson](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/01_Analytical_solutions/2_Isotropic_elastic_medium.ipynb) the momentum conservation and stress-strain relation for the isotropic 3D elastic medium in Cartesian coordinates can be written in **stress-velocity formulation** as:\n",
    "\n",
    "\\begin{align}\n",
    "\\rho \\dot{v}_x &= \\frac{\\partial \\sigma_{xx}}{\\partial x} + \\frac{\\partial \\sigma_{xy}}{\\partial y} + \\frac{\\partial \\sigma_{xz}}{\\partial z} + f_x\\; \\hspace{1.5cm} \\dot{\\sigma}_{xx} = \\lambda \\dot{\\Delta} + 2 \\mu \\dot{\\epsilon}_{xx}\\nonumber \\\\\n",
    "\\rho \\dot{v}_y &= \\frac{\\partial \\sigma_{yx}}{\\partial x} + \\frac{\\partial \\sigma_{yy}}{\\partial y} + \\frac{\\partial \\sigma_{yz}}{\\partial z} + f_y\\; \\hspace{1.5cm}\n",
    "\\dot{\\sigma}_{yy} = \\lambda \\dot{\\Delta} + 2 \\mu \\dot{\\epsilon}_{yy}\n",
    "\\nonumber \\\\\n",
    "\\rho \\dot{v}_z &= \\frac{\\partial \\sigma_{zx}}{\\partial x} + \\frac{\\partial \\sigma_{zy}}{\\partial y} + \\frac{\\partial \\sigma_{zz}}{\\partial z} + f_z\\; \\hspace{1.5cm} \n",
    "\\dot{\\sigma}_{zz} = \\lambda \\dot{\\Delta} + 2 \\mu \\dot{\\epsilon}_{zz}\\nonumber \\\\\n",
    "\\dot{\\sigma}_{xy} &= 2 \\mu \\dot{\\epsilon}_{xy}\\; \\hspace{1.45cm} \\dot{\\sigma}_{xz} = 2 \\mu \\dot{\\epsilon}_{xz}\\; \\hspace{1.45cm} \\dot{\\sigma}_{yz} = 2 \\mu \\dot{\\epsilon}_{yz}\\nonumber \\\\\n",
    "\\end{align}\n",
    "\n",
    "with \n",
    "\n",
    "\\begin{align}\n",
    "v_{i} &\\; \\hspace{0.1cm} \\text{particle velocity vector}\\; \\hspace{1.5cm} \\epsilon_{ij} \\; \\hspace{0.1cm} \\text{strain tensor}\\notag\\\\\n",
    "f_{i} &\\; \\hspace{0.1cm} \\text{source vector}\\; \\hspace{2.9cm} \\lambda,\\; \\mu \\; \\hspace{0.1cm} \\text{Lamé parameters}\\notag\\\\\n",
    "\\sigma_{ij} &\\; \\hspace{0.1cm} \\text{stress tensor}\\; \\hspace{3.6cm} \\rho \\; \\hspace{0.1cm} \\text{density}\\notag\\\\\n",
    "\\dot{\\Delta} & = \\dot{\\epsilon}_{xx}+\\dot{\\epsilon}_{yy}+\\dot{\\epsilon}_{zz}\\notag \\; \\hspace{2.5cm} \\dot{\\epsilon}_{ij} = \\frac{1}{2}\\biggl(\\frac{\\partial v_i}{\\partial x_j} + \\frac{\\partial v_j}{\\partial x_i}\\biggr)\n",
    "\\end{align}\n",
    "\n",
    "Introducing a spherical coordinate system:\n",
    "\n",
    "<img src=\"images/3D_spherical.png\" width=\"30%\">\n",
    "\n",
    "where $\\varphi$ denotes the azimuth angle, $\\theta$ the polar angle, $r$ the radius, we can transform the partial differential equations from the cartesian to spherical coordinate system. The easiest approach is to write the governing equations in vector-tensor notation:\n",
    "\n",
    "\\begin{align}\n",
    "\\rho \\frac{\\partial \\mathbf{v}}{\\partial t} &= \\mathbf{\\nabla} \\cdot \\mathbf{\\sigma} + \\mathbf{f} \\nonumber \\\\\n",
    "\\frac{\\partial \\mathbf \\sigma}{\\partial t} &= \\mathbf{c} : \\frac{\\partial \\mathbf \\epsilon}{\\partial t} \\nonumber \\\\\n",
    "\\frac{\\partial \\mathbf \\epsilon}{\\partial t} &= \\frac{1}{2}\\biggl(\\mathbf{\\nabla\\; v} + (\\mathbf{\\nabla\\; v})^T\\biggr) \\nonumber\n",
    "\\end{align}\n",
    "\n",
    "and evaluate the divergence and gradient operators in spherical coordinates. For example, the divergence operator applied to the 2nd order stress tensor $\\sigma_{ij}$ in spherical coordinates is defined [here](https://de.wikipedia.org/wiki/Formelsammlung_Tensoranalysis#Divergenz_in_verschiedenen_Koordinatensystemen). You can find a similar definition for the gradient applied to the vector $v_i$. \n",
    "\n",
    "##### Bonus exercise\n",
    "\n",
    "Use the definition of a gradient applied to the vector $v_i$ in spherical coordinates to compute the time derivative of the strain tensor $\\dot{\\epsilon}$.\n",
    "\n",
    "The resulting equations of motion and stress-strain relations in spherical coordinates are ...\n",
    "\n",
    "\\begin{align}\n",
    "\\rho \\dot{v}_r &= \\frac{\\partial \\sigma_{rr}}{\\partial r} + \\frac{1}{r}\\frac{\\partial \\sigma_{r\\theta}}{\\partial \\theta} + \\frac{1}{r \\sin{\\theta}} \\frac{\\partial \\sigma_{r\\varphi}}{\\partial \\varphi} + \\frac{1}{r}(2 \\sigma_{rr} - \\sigma_{\\theta\\theta} - \\sigma_{\\varphi\\varphi} + \\sigma_{r \\theta}\\; cot \\theta) + f_r\\nonumber \\\\\n",
    "\\rho \\dot{v}_\\theta &= \\frac{\\partial \\sigma_{r\\theta}}{\\partial r} + \\frac{1}{r}\\frac{\\partial \\sigma_{\\theta\\theta}}{\\partial \\theta} + \\frac{1}{r \\sin{\\theta}} \\frac{\\partial \\sigma_{\\theta\\varphi}}{\\partial \\varphi} + \\frac{1}{r}\\{(\\sigma_{\\theta\\theta} - \\sigma_{\\varphi\\varphi})\\;cot \\theta + 3 \\sigma_{r\\theta}\\} + f_\\theta\\nonumber \\\\\n",
    "\\rho \\dot{v}_\\varphi &= \\frac{\\partial \\sigma_{r\\varphi}}{\\partial r} + \\frac{1}{r}\\frac{\\partial \\sigma_{\\theta\\varphi}}{\\partial \\theta} + \\frac{1}{r \\sin{\\theta}} \\frac{\\partial \\sigma_{\\varphi\\varphi}}{\\partial \\varphi} + \\frac{1}{r}(3 \\sigma_{r\\varphi} + 2\\sigma_{\\theta\\varphi}\\; cot \\theta) + f_\\varphi\\nonumber \\\\\n",
    "& \\notag\\\\\n",
    "\\dot{\\sigma}_{rr} &= \\lambda \\dot{\\Xi} + 2 \\mu \\dot{\\epsilon}_{rr}\\; \\hspace{1.45cm} \\dot{\\sigma}_{\\theta\\theta} = \\lambda \\dot{\\Xi} + 2 \\mu \\dot{\\epsilon}_{\\theta\\theta}\\; \\hspace{1.45cm} \\dot{\\sigma}_{\\varphi\\varphi} = \\lambda \\dot{\\Xi} + 2 \\mu \\dot{\\epsilon}_{\\varphi\\varphi}\\nonumber \\\\\n",
    "\\dot{\\sigma}_{r\\theta} &= 2 \\mu \\dot{\\epsilon}_{r\\theta}\\; \\hspace{2.42cm} \\dot{\\sigma}_{\\theta\\varphi} = 2 \\mu \\dot{\\epsilon}_{\\theta\\varphi}\\; \\hspace{2.55cm} \\dot{\\sigma}_{r\\varphi} = 2 \\mu \\dot{\\epsilon}_{r\\varphi}\\nonumber \\\\\n",
    "& \\notag \\\\\n",
    "\\dot{\\epsilon}_{rr} &= \\frac{\\partial v_r}{\\partial r}\\; \\hspace{1.45cm} \n",
    "\\dot{\\epsilon}_{\\theta\\theta} = \\frac{1}{r}\\frac{\\partial v_\\theta}{\\partial \\theta} + \\frac{1}{r} v_r\\; \\hspace{1.45cm}\n",
    "\\dot{\\epsilon}_{\\varphi\\varphi} = \\frac{1}{r\\; sin\\theta}\\frac{\\partial v_\\varphi}{\\partial \\varphi} + \\frac{1}{r} v_r + \\frac{cot \\theta}{r} v_\\theta\\; \\hspace{1.45cm}\\nonumber \\\\\n",
    "\\dot{\\epsilon}_{r\\theta} &= \\frac{1}{2}\\biggl(\\frac{1}{r} \\frac{\\partial v_r}{\\partial \\theta} + \\frac{\\partial v_\\theta}{\\partial r} - \\frac{1}{r} v_\\theta\\biggr)\\; \\hspace{2.6cm} \\dot{\\epsilon}_{r\\varphi} = \\frac{1}{2}\\biggl(\\frac{1}{r\\; sin \\theta} \\frac{\\partial v_r}{\\partial \\varphi} + \\frac{\\partial v_\\varphi}{\\partial r} - \\frac{1}{r} v_\\varphi\\biggr)\n",
    "\\nonumber \\\\\n",
    "\\dot{\\epsilon}_{\\theta\\varphi} &= \\frac{1}{2}\\biggl(\\frac{1}{r\\; sin\\theta}\\frac{\\partial v_\\theta}{\\partial \\varphi} + \\frac{1}{r} \\frac{\\partial v_\\varphi}{\\partial \\theta} - \\frac{cot \\theta}{r} v_{\\varphi} \\biggr)\\nonumber \n",
    "\\end{align}\n",
    "\n",
    "with \n",
    "\n",
    "$\\Xi = \\dot{\\epsilon}_{rr} + \\dot{\\epsilon}_{\\theta\\theta} + \\dot{\\epsilon}_{\\varphi\\varphi}$\n",
    "\n",
    "Puuh, that problem looks quite time consuming in terms of writing a FD code as well as computation time. Furthermore, there exist singularities at the center of the planet  ($r=0$) due to the $1/r$ terms and at the poles ($\\theta=0,\\; \\pi$) related to the $cot \\theta$- and $1/sin(\\theta)$-terms. In order to model global seismic wave propagation on your gaming laptop in a reasonable time, we have to simplify the problem ...\n",
    "\n",
    "## Axisymmetric SH problem in spherical coordinates\n",
    "\n",
    "Similar to the SH problem in Cartesian coordinates, we can derive the SH problem in spherical coordinates by assuming that $\\lambda=0$ and the only non-zero particle velocity component is $v_\\varphi \\neq 0$, while $v_\\theta = v_r = 0$. This simplifies the 3D problem to:\n",
    "\n",
    "\\begin{align}\n",
    "\\rho \\dot{v}_\\varphi &= \\frac{\\partial \\sigma_{r\\varphi}}{\\partial r} + \\frac{1}{r}\\frac{\\partial \\sigma_{\\theta\\varphi}}{\\partial \\theta} + \\frac{1}{r \\sin{\\theta}} \\frac{\\partial \\sigma_{\\varphi\\varphi}}{\\partial \\varphi} + \\frac{1}{r}(3 \\sigma_{r\\varphi} + 2\\sigma_{\\theta\\varphi}\\; cot \\theta) + f_\\varphi\\nonumber \\\\\n",
    "& \\notag\\\\\n",
    "\\dot{\\sigma}_{\\varphi\\varphi} &= 2 \\mu \\dot{\\epsilon}_{\\varphi\\varphi}\\; \\hspace{2.55cm}\n",
    "\\dot{\\sigma}_{r\\varphi} = 2 \\mu \\dot{\\epsilon}_{r\\varphi}\\; \\hspace{3.65cm} \\dot{\\sigma}_{\\theta\\varphi} = 2 \\mu \\dot{\\epsilon}_{\\theta\\varphi}\\nonumber \\\\\n",
    "& \\notag \\\\\n",
    "\\dot{\\epsilon}_{\\varphi\\varphi} &= \\frac{1}{r\\; sin\\theta}\\frac{\\partial v_\\varphi}{\\partial \\varphi}\\; \\hspace{1.6cm} \\dot{\\epsilon}_{r\\varphi} = \\frac{1}{2}\\biggl(\\frac{\\partial v_\\varphi}{\\partial r} - \\frac{1}{r} v_\\varphi\\biggr)\\; \\hspace{1.45cm} \\dot{\\epsilon}_{\\theta\\varphi} = \\frac{1}{2}\\biggl(\\frac{1}{r} \\frac{\\partial v_\\varphi}{\\partial \\theta} - \\frac{cot \\theta}{r} v_{\\varphi} \\biggr)\\nonumber \n",
    "\\end{align}\n",
    "\n",
    "The problem can be further simplified by assuming that the wavefields and model parameters are invariant in $\\varphi$-direction, which means that all spatial derivatives $\\frac{\\partial}{\\partial \\varphi} = 0$, leading to the **axisymmetric spherical SH problem**\n",
    "\n",
    "\\begin{align}\n",
    "\\rho \\dot{v}_\\varphi &= \\frac{\\partial \\sigma_{r\\varphi}}{\\partial r} + \\frac{1}{r}\\frac{\\partial \\sigma_{\\theta\\varphi}}{\\partial \\theta} + \\frac{1}{r}(3 \\sigma_{r\\varphi} + 2\\sigma_{\\theta\\varphi}\\; cot \\theta) + f_\\varphi\\nonumber \\\\\n",
    "& \\notag\\\\\n",
    "\\dot{\\sigma}_{r\\varphi} &= 2 \\mu \\dot{\\epsilon}_{r\\varphi}\\; \\hspace{3.65cm} \\dot{\\sigma}_{\\theta\\varphi} = 2 \\mu \\dot{\\epsilon}_{\\theta\\varphi} \\\\\n",
    "& \\notag \\\\\n",
    "\\dot{\\epsilon}_{r\\varphi} &= \\frac{1}{2}\\biggl(\\frac{\\partial v_\\varphi}{\\partial r} - \\frac{1}{r} v_\\varphi\\biggr)\\; \\hspace{1.45cm} \\dot{\\epsilon}_{\\theta\\varphi} = \\frac{1}{2}\\biggl(\\frac{1}{r} \\frac{\\partial v_\\varphi}{\\partial \\theta} - \\frac{cot \\theta}{r} v_{\\varphi} \\biggr)\\nonumber \n",
    "\\end{align}\n",
    "\n",
    "Compared to the Cartesian 2D SH problem, the axisymmetric SH problem describes wave propagation with correct 3D geometrical spreading with the limitation that the model has to be axisymmetric. Another disadvantage is that we can only model ring-sources, if the source is not placed on the symmetry axis, which is actually not possible due to the $cot \\theta$ term.\n",
    "\n",
    "## Finite-Difference solution of the axisymmetric SH problem\n",
    "\n",
    "We can apply a similar staggered grid finite-difference approach for the axisymmetric SH problem eqs. (1) as we have used for the [Cartesian SH problem](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/06_2D_SH_Love_wave_modelling/1_2D_SH_FD_staggered.ipynb)\n",
    "\n",
    "<img src=\"images/SG_SH-sph.png\" width=\"50%\">\n",
    "\n",
    "The correct FD operators, material parameters and coordinates on the staggered grid are estimated in two steps:\n",
    "\n",
    "1. Where are the material parameters and wavefields on the LHS of the equation located on the staggered grid?\n",
    "\n",
    "2. All derivatives on the RHS have to be centered on the grid point estimated in step 1. Depending on the equation you also have to define material parameters, wavefields and coordinates on the correct positions.\n",
    "\n",
    "Applying these rules, we can estimate the following 2nd order spatial/temporal staggered grid FD scheme with explicit Leap frog time stepping for the momentum equation:\n",
    "\n",
    "\\begin{align}\n",
    "\\rho(i,j) \\frac{v_{\\phi}(i,j,n+1/2) - v_{\\phi}(i,j,n-1/2)}{dt} &= \\partial_r \\sigma_{r\\varphi} + \\frac{1}{r} \\partial_\\theta \\sigma_{\\theta\\varphi} + \\frac{1}{r}(3 \\sigma_{r\\varphi} + 2\\sigma_{\\theta\\varphi}\\; cot \\theta) + f_\\varphi\\notag\n",
    "\\end{align}\n",
    "\n",
    "with the spatial approximations:\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_r \\sigma_{r\\varphi} &\\approx \\frac{\\sigma_{r\\varphi}(i,j+1/2,n)-\\sigma_{r\\varphi}(i,j-1/2,n)}{dr}\\notag\\\\\n",
    "\\frac{1}{r} \\partial_\\theta \\sigma_{\\theta\\varphi} &\\approx \\frac{1}{r(i,j)}\\frac{\\sigma_{\\theta\\varphi}(i+1/2,j,n)-\\sigma_{\\theta\\varphi}(i-1/2,j,n)}{d\\theta}\\notag\\\\\n",
    "\\frac{3}{r} \\sigma_{r\\varphi} & \\approx \\frac{3}{r(i,j)} \\frac{\\sigma_{r\\varphi}(i,j+1/2,n)+\\sigma_{r\\varphi}(i,j-1/2,n)}{2}\\notag\\\\\n",
    "\\frac{2}{r} \\sigma_{\\theta\\varphi} cot\\; \\theta & \\approx \\frac{2}{r(i,j)}\\frac{\\sigma_{\\theta\\varphi}(i+1/2,j,n)+\\sigma_{\\theta\\varphi}(i-1/2,j,n)}{2} cot\\; \\theta(i,j)\\notag\\\\\n",
    "f_\\varphi &\\approx f_\\varphi(i,j) \\notag\n",
    "\\end{align}\n",
    "\n",
    "and stress-strain relation:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\sigma_{r\\varphi}(i,j+1/2,n+1)-\\sigma_{r\\varphi}(i,j+1/2,n)}{dt} &= \\mu \\biggl(\\partial_r v_\\varphi - \\frac{1}{r} v_\\varphi\\biggr) \\notag\\\\\n",
    "\\frac{\\sigma_{\\theta\\varphi}(i+1/2,j,n+1)-\\sigma_{\\theta\\varphi}(i+1/2,j,n)}{dt} &= \\mu  \\biggl(\\frac{1}{r} \\partial_\\theta v_\\varphi - \\frac{cot \\theta}{r} v_{\\varphi} \\biggr)\\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "with\n",
    "\n",
    "\\begin{align}\n",
    "\\mu \\partial_r v_\\varphi &\\approx \\mu_r(i,j+1/2) \\frac{v_{\\varphi}(i,j+1,n+1/2)-v_{\\varphi}(i,j,n+1/2)}{dr}\\notag\\\\\n",
    "\\mu \\frac{1}{r} v_\\varphi &\\approx \\mu_r(i,j+1/2) \\frac{1}{r(i,j+1/2)} \\frac{v_{\\varphi}(i,j+1,n+1/2)+v_{\\varphi}(i,j,n+1/2)}{2}\\notag\\\\\n",
    "\\mu \\frac{1}{r}\\partial_\\theta v_\\varphi &\\approx \\mu_\\theta(i+1/2,j) \\frac{1}{r(i+1/2,j)}\\frac{v_{\\varphi}(i+1,j,n+1/2)-v_{\\varphi}(i,j,n+1/2)}{d\\theta}\\notag\\\\\n",
    "\\mu \\frac{cot\\; \\theta}{r} v_\\varphi &\\approx \\mu_\\theta(i+1/2,j) \\frac{cot\\; \\theta(i+1/2,j)}{r(i+1/2,j)}\\frac{v_{\\varphi}(i+1,j,n+1/2)+v_{\\varphi}(i,j,n+1/2)}{2}\\notag\n",
    "\\end{align}\n",
    "\n",
    "with the harmonically averaged shear moduli:\n",
    "\n",
    "\\begin{align}\n",
    "\\mu_r(i,j+1/2) &= 2 [1/\\mu(i,j)+1/\\mu(i,j+1)]^{}-1 \\notag\\\\\n",
    "\\mu_\\theta(i+1/2,j) &= 2 [1/\\mu(i,j)+1/\\mu(i+1,j)]^{}-1 \\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "### OpenSource codes to model seismic wave propagation in spherical coordinates\n",
    "\n",
    "In the next Jupyter notebook we will implement the FD solution for the axisymmetric SH problem. However, there exist multiple OpenSource codes, who deal with this problem:\n",
    "\n",
    "- [**SHaxi**](http://www.quest-itn.org/library/software/shaxi) by Gunnar Jahnke, Michael Thorne and Heiner Igel is basically a Fortran 90 implementation of the FD scheme described above. Furthermore, it is parallelized by domain decomposition allowing it to run on high performance parallel computers for global seismic wave propagation\n",
    "\n",
    "- [**AxiSEM**](http://seis.earth.ox.ac.uk/axisem/) by Martin van Driel, Simon Stähler, Tarje Nissen-Meyer, Stefanie Hempel and Alexandre Fournier solves the more general axisymmetric (an)isotropic (visco)-elastic problem using the Spectral Element method. The code is parallelized by domain decomposition using MPI\n",
    "\n",
    "- [**SPECFEM3D Globe**](https://specfem.org/) originally developed by Dimitri Komatitsch and Jean-Pierre Vilotte, now supported by a broad community, solves the general 3D (an)isotropic (visco)-elastic in spherical coordinates using Spectral Elements for modelling global wave propagation. Beside MPI parallelization it also supports CUDA to run SPECFEM3D Globe on GPUs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We learned:\n",
    "\n",
    "* 3D equations of motion for an isotropic elastic medium in spherical coordinates\n",
    "* Introducing the 2D axisymmetric SH approximation by defining $v_\\varphi$ as the only non-zero particle velocity component and assuming invariance of wavefields and material parameters in $\\varphi$-direction\n",
    "* Staggered grid FD solution of the axisymmetric SH problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "1. Jahnke, G. (2009) \"[Methods for Seismic Wave Propagation on Local and Global Scales with Finite Differences](https://edoc.ub.uni-muenchen.de/11235/)\", PhD thesis, Ludwig-Maximilians University Munich."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
